Models for the evolution of free granular 

surfaces 

R. Mulet* and H. Herrmann 
ICAl, University of Stuttgart, 
Pfaffenwaldring 27, 70596 Stuttgart, 
Germany 

February 1, 2008 



Abstract 

We introduce two sets of continuum equations to describe granular 
flow on a free surface and study their properties. The equations de- 
rived from a microscopic picture that includes jumps and a mobility 
threshold, account for ripple and crater formation. 
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1 Introduction 

One of the apparently simple problems for which there is no consensus yet on 
which are the correct continuum equations of motion is the flux of grains on 
a granular surface (see for example ref. Jl], ^ §]). This flux can be driven by 
a fluid, like air, by gravity or by an initial impact. It modifies the granular 
surface itself either by deposition or by erosion. We will consider in this paper 
a gravity driven flow with variable initial energy and consider deposition only. 
In addition we want to be in the limit of small flux. One realization is the 
slow buildup of a pile from a point source. 
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A discrete model was proposed in ref. in which particles jump down 
stairs and deposit if their energy is below a certain material dependent thresh- 
old. From this model a set of continuum equations was derived in ref. |^ and 
applied to the problem of impact craters. This continuum description, how- 
ever, made an assumption on the time dependence of the slope that only 
contained the threshold in an indirect way, namely by inserting the known 
steady state result. Also all the jumps were of length unity even if their 
kinetic energy was different. In this paper we want to deal with these defects 
of the equations in ref. P| by inserting the two mentioned effects indepen- 
dently into the equations and investigating the consequences analytically and 
numerically. 

2 Modelling jumps of grains 

We will introduce an extension to the model proposed in ref. There, 
particles of mass m are added at the top of a pile, and they move, under 
the action of gravity g, until their kinetic energy e falls below a certain 
threshold U (related to the friction between grains) where they then stop. 
At each collision the kinetic energy decreases by a factor proportional to the 
restitution coefficient of the material r. Mathematically, it reads: 

e{x + S) = (e(x) + mg{h{x) — h{x + S)))r (1) 

It was assumed in ref. that the traveled distance after each particle 
jump is constant. Here, trying to model a more realistic situation, we assume 
that this distance is proportional to the kinetic energy of the particle. So, 
we propose to modify equation (|1|) to: 

e(x + qe) = {e{x) + mg{h{x) — h{x + qe)))r (2) 

where g is a proportionality coefficient (in units of ^). Then, writing eq. 
(H) in differential form we obtain: 

^ = -(r - 1) + r7(x) (3) 
ax q 

where 7(x) = —mgdh/dx. 
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To this equation one adds in ref . P an evolution equation for the slope of 
the form: 



g = r(f/-e(x)) (4) 

where F controls the rate of relaxation of the slope. Boundary conditions were 
imposed trying to reproduce real experimental situations of heap formation: 
e(0) = Co a value which is proportional to the height from which the grains 
fall on the pile, and 7(0) = 0. 

The equations and (H), constitute a pair of linear differential equations 
with constant coefficients. From them, the angle of a pile can be easily 
calculated: 

7(x) = (1 — cos(vTrx)) — vTrfco — U) sm( VTrx) (5) 

rq 

and then, integrating equation (Bf) the shape of the pile is: 



1 — r 1 — r 1 / — 1 / — 

h(x) = h(0) x-\ =sin(vFra;) (co — f/) cos(v Frx)) (6) 

rqmg rqmg y Fr rmg 

This solution has two remarkable properties as shown in Figure 1. First 
of all the oscillatory behaviour superimposed to the usual linear decay of the 
pile. These oscillations arise from the existence of a length scale F which 
controles the rate of the relaxation of the slope while independently the 
particles jump a distance (ge) determined only by their kinetic energy. This 
phenomenon is similar to the formation of ripples due to saltation , and it 
could be the first theoretical prediction of "gravity" induced ripples. 

Another important feature of eq. (6) is that it can reproduce a crater 
formation on the top of pile. In fact, in the hmit vTrx << 1, equation (|^) 
transforms into: 

mgh{x) = mgh{Q) - 1 (e„ - f/) (1 - ^-x") (7) 

showing the existence of a parabolic-like crater whose depth linearly increases 
with the initial kinetic energy of the particles. 
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3 Modelling the mobility threshold 

In this model, we add to the equation for the energies from ref. 

de r — 1 , , 

- = ^e + r7(x) (8) 

a new equation for the angles based on the assumption that the equilibrium 
local slope of the pile is constant, i.e. 

0(x + 6) = (9) 

where is the slope of the pile. 

The evolution of the slope proceeds in the following way: if a particle 
arrives at x, it has two possibilities (as described for a similar model in ref. 
1^]), if e > ?7, it moves to x + 6, then: 

(j){x + 6) = (j){x) - 1 (10) 

or, it stops when e < U and 

(f){x + 6) = (f){x) + 1 (11) 
Equations (|10D and may be written in the following differential form: 

^ = ]'9n{U - e) (12) 

We next solve the coupled differential equations (||) and (|12D numerically. 
It is more convenient to write both equations in the same form. So dividing 
eq. (IP by mg, and noting that 7 = mgcj), we obtain the following equations: 

dz r — 1 r ^ 

_ = -sgn{zu - z) (14) 

where mgz = e and mgzu = U. 

Then, once (pix) is known the height of the pile can be calculated through 
integration of (p and obtained numerically. 



To numerically solve the equations (0) and (|Tj) the function sgn[x 



was replaced by the continous function tanh{ax). Profiles calculated using 
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different values of a are shown in Figure 2. From the Figure we conclude 
that if a > 10^ all the numerical results are equivalent. To be on the secure 
side we used in the rest of the paper a = 10^. 

Figure 3 represents typical sandpile profiles obtained for different values 
of the parameters. The bottom curve (where the crater is better defined) 
represents a system where the grains have higher initial energy. 

In figure 4 is plotted the equilibrium angle of the pile (j)eq as a function of 
U{1 — r)/r for different values of Cq. As previously calculated in ref. 0], the 
angle of the pile is equal to f/(l — r)/r showing the equivalence between our 
approach and that proposed in [|]. 

The depth Ah and the width Ax of the craters, should only depend on 
the three parameters involved in the model r, Co and U. In fact, in figures 
5 and 6, Ah and Ax obtained using different sets of parameters (r, Cq, U) 
collapse on the same curve following the scaling relations: 

Ah = C^-irf{C^-l)-^r) (15) 

and 

Ax=(^-l)V(^-irr) (16) 

where (3 = 2.0, a = 1.0, u = 4.0 and f{x) and g{x) are scaling functions with 
the following properties f{x) ~ 1 if x << 1 and f{x) ~ if x >> 1, idem for 
g. This indicates that around Cq = U and r = we have critical behaviour 
with simple integer exponents. 

Physically equations (|T5[) and (|1^) imply that in granular systems with 
negligible restitution coefficent (r ~ 0) the depth of the crater increases 
linearly with Co/U — 1 and with 1/r, in good agreement with eq.|^ while the 
width of the craters increases parabolically [jS = 2.0). Therefore systems 
with small Cq, i.e. (cq ~ f/) or r ~ 1 do not develop a crater. 

4 Conclusions 

We have presented two sets of equations of motion for the limit of dilute 
granular surface flow. The first set included a jump length proportional to 
the energy of the grains. This gave ripple formation on the critical slope, a 
phenomenon which has not yet been observed experimentally. The second 
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set was a cleaner way to include the deposition threshold than done in ref. 
[0]. The results confirm that the previously used equation (introduced in ref. 
IQ]) gave qualitatively the right answer for the shape of the pile which means 
that the picture of a characteristic length for relaxation to the equilibrium 
angle is correct. 
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Figure Captions 



Figure 1 Graphical representation of equation 6. We chose = 

0.1, VrV = 1 and = 0.1. 

Figure 2 Shape of the pile calculated by integrating the solution of 
equations (13) and (14) for a = 1 (upper curve) and a = 10 and 100 (bottom 
curve). 

Figure 3 Shape of the pile for r = 0.03, zu = 0.05, and (from top to 
bottom) z = 0.05,0.10 and 15. The profiles are shifted in the z axis to 
distinguish them better. 

Figure 4 Equilibrium angle of the pile (peg as a function of f/(l — r)/r 
for different values of eo,r and U. 

Figure 5 Data collapse of the depth of the craters using eq: (p!5|) .Parameters: 
(r,eo,f/) = A (0.01, 0.15, 0.05); □ (0.01, 0.25, 0.05); (0.01, 0.25, 0.10); 
• (0.01, 0.25, 0.15); x (0.03, 0.20, 0.10); + (0.03, 0.25, 0.15); □ (0.05, 0.30, 
0.10). 

Figure 6 Data collapse of the depth of the crater using eq: (|16]). Param- 
eters: (r,eo,f/) = A (0.01, 0.15, 0.05); □ (0.01, 0.25, 0.05); (0.01, 0.25, 
0.10); • (0.01, 0.25, 0.15); x (0.03, 0.20, 0.10); + (0.03, 0.25, 0.15); □ (0.05, 
0.30, 0.10). 
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